import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score, mean_squared_error
import math

# 设置随机种子以确保结果可复现
np.random.seed(42)

# 生成参考浓度（真值）数据点
true_concentrations = np.array([
    0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100,
    0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100
])

# 构建模型预测结果（带一些随机误差）
predicted_concentrations = true_concentrations * (1 + 0.1 * np.random.normal(size=len(true_concentrations)))

# 区分建模集和独立预测集
modeling_set_indices = np.array([True] * 11 + [False] * 11)
prediction_set_indices = ~modeling_set_indices

# 计算预测集性能指标
rmsep = math.sqrt(mean_squared_error(
    true_concentrations[prediction_set_indices],
    predicted_concentrations[prediction_set_indices]
))
r_squared = r2_score(
    true_concentrations[prediction_set_indices],
    predicted_concentrations[prediction_set_indices]
)

# 创建散点图
plt.figure(figsize=(10, 6))

# 绘制建模集数据点（蓝色圆圈）
plt.scatter(
    np.log10(true_concentrations[modeling_set_indices] + 1),
    np.log10(predicted_concentrations[modeling_set_indices] + 1),
    color='blue',
    marker='o',
    label='Modeling Set',
    alpha=0.7
)

# 绘制预测集数据点（红色三角）
plt.scatter(
    np.log10(true_concentrations[prediction_set_indices] + 1),
    np.log10(predicted_concentrations[prediction_set_indices] + 1),
    color='red',
    marker='^',
    label='Prediction Set',
    alpha=0.7
)

# 添加理想拟合线（Y=X线）
lims = [
    np.min(np.log10(true_concentrations + 1)) - 0.1,
    np.max(np.log10(true_concentrations + 1)) + 0.1
]
plt.plot(lims, lims, 'k--', alpha=0.7, label='Y=X Line')

# 标注模型性能指标
plt.text(
    lims[0] + 0.1,
    lims[1] - 0.1,
    f'R² = {r_squared:.3f}\nRMSEP = {rmsep:.2f} ppm',
    fontsize=10,
    bbox=dict(facecolor='white', alpha=0.7)
)

# 设置坐标轴标签
plt.xlabel('Reference Concentration (ppm) [log10(Conc+1)]')
plt.ylabel('Predicted Concentration (ppm) [log10(Conc+1)]')
plt.title('PLSR Prediction Results')

# 添加图例
plt.legend()

# 设置对数坐标轴刻度标签
plt.xticks(ticks=np.log10(np.array([0.1, 1, 10, 100]) + 1),
           labels=['0.1', '1', '10', '100'])
plt.yticks(ticks=np.log10(np.array([0.1, 1, 10, 100]) + 1),
           labels=['0.1', '1', '10', '100'])

# 显示散点图
plt.show()